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Abstract 

We give a brief discussion of the recently developed Constrained-Path Monte Carlo 
Method. This method is a quantum Monte Carlo technique that eliminates the 
fermion sign problem plaguing simulations of systems of interacting electrons. The 
elimination is accomplished by trading an exact procedure for an approximate one 
that has been demonstrated to give very accurate estimates of energies and many- 
body correlation functions. We also give a short review of its applications, a dis- 
cussion of several strategies for parallelizing it, and some speculation of its future 
extensions. 



1 Introduction 



We will give a brief discussion of the recently developed Constrained-Path 
Monte Carlo (CPMC) Method [1,2]. This ground-state (T = 0) method is a 
quantum Monte Carlo (QMC) technique that eliminates the infamous fermion 
sign problem plaguing simulations of systems of interacting electrons. The 
fermion sign problem causes the variance of measured quantities to increase 
exponentially with increasing system size and decreasing temperature. Rapidly 
the sign problem destroys one's ability to compute with acceptable accuracy. 
With the CPMC method one has simulated system sizes not possible with the 
standard method [3]. In particular, the lattice size dependence of many-body 
superconducting pairing correlations functions were simulated for some of the 
largest lattice sizes studied to date. 

In CPMC method the elimination of the sign problem is accomplished by 
trading an exact procedure for an approximate one that has been demon- 
strated to give very accurate estimates of energies and various many-body 
correlation functions. The exact procedure determines the lowest eigenvalue 
and eigenvector of the Hamiltonian by projecting them from a trial state. This 
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procedure is easily converted to a branched random walk. Because of the sign 
problem, the random walkers carry a positive and negative weight, unfortu- 
nately in proportions such that the average weight (sign) becomes zero as the 
system size increases. The constrained path method is a particular way to 
break the symmetry in the sign of the walkers and produce only positively 
weighted walkers by eliminating those with a negative overlap with a certain 
constraining state. The procedure bears a similarly to the fixed-node Monte 
Carlo method [4] that has been successfully used for several decades in simula- 
tions of interacting electrons systems defined in the continuum of configuration 
space. The CPMC method however operates in the manifold of single-particle 
states (Slater determinants) defined in Fock space and hence represents quite 
a different and novel Monte Carlo algorithm. 

In the next section we will give a brief discussion of the method. Then in the 
following sections will discuss the various models to which the method has 
been applied, highlighting significant results. After this we will discuss several 
strategies for parallelizing the method. At first glance it would seem as if we 
could simply exploit the natural parallelization enjoyed by most Monte Carlo 
methods [5] . In fact we need to do a bit more. In the closing section we make 
some speculations on future applications and extensions of the CPMC method 
and the potential changes in parallelization procedures. 



2 The Constrained-Path Monte Carlo Method 



Our numerical method is extensively described and benchmarked elsewhere 
[1,2]. Here we only discuss its basic strategy and approximation. In the CPMC 
method, the ground-state wave function 1^*0) is projected from a known initial 
wave function |\E'r) by a branching random walk in an over-complete space 
of Slater determinants In such a space, we can write l^'o) = S</)X(0)|0)- 
The random walk produces an ensemble of called random walkers, which 
represent |\t'o) in the sense that their distribution is a Monte Carlo sampling 
of x(0), that is, a sampling of the ground-state wave function. 

To completely specify the ground-state wave function for a system of inter- 
acting electrons, only determinants satisfying (\E'o|0) > are needed because 
1^0) resides in either of two degenerate halves of the Slater determinant space, 
separated by a nodal surface N that is defined by {'^o\4>) = 0. The degeneracy 
is a consequence of both \ipo) and — |^o) satisfying Schrodinger's equation. The 
sign problem occurs because walkers can cross N as their orbitals evolve con- 
tinuously in the random walk. Asymptotically they populate the two halves 
equally, leading to an ensemble that has zero overlap with |\E'o)- If N were 
known, we would simply constrain the random walk to one half of the space 
and obtain an exact solution of Schrodinger's equation. In the constrained- 
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path QMC method, without a priori knowledge of N, we use a trial wave 
function |^7-) and require > 0. This is what is called the constrained- 

path approximation. 

The quality of the calculation clearly depends on the trial wave function |^'t)- 
Since the constraint only involves the overall sign of its overlap with any deter- 
minant 10), it seems reasonable to expect the results to show some insensitivity 
to |\£'r)- Through extensive benchmarking on the Hubbard model, it has been 
found that simple choices of this function can give very good results [1,2]. 

Besides as starting point and as a condition constraining a random walker, we 
also use \^t) as an importance function. To reduce variance, we use (^t|0) 
to bias the random walk into those parts of Slater determinant space that 
have a large overlap with the trial state. For all three uses of |^t)) it clearly 
is advantageous to have |\E'r) approximate |\E'o) as closely as possible. Only in 
the constraining of the path does |\l'r) 7^ |^o) generate an approximation. 

Almost all the calculations reported here are done for square lattices with 
periodic boundary conditions. Mostly, we study closed shell cases, for which 
the corresponding free-electron wave function is non-degenerate and is transla- 
tionally invariant. In these cases, the free-electron wave function, represented 
by a single Slater determinant, is used as the trial wave function \iPt)- The 
use of an unrestricted Hartree-Fock wave function as I'^t) generally produces 
no significant improvement in the results. 

Wc remark the the CPMC method has been extended to use generalized 
Hartree-Fock wave functions, of which the most famous example is the BCS 
wave function. The trick for doing this is described elsewhere [6]. 

One does simulations because exact solutions are generally unavailable, and 
approximate solutions, like those from Hartree-Fock approximations, are of- 
ten poor. The objective of the CPMC method and other simulation methods 
is to go beyond Hartree Fock approximation. By expressing the wave func- 
tion as a hnear combination of Slater determinant, the CPMC method is a 
type of stochastic configuration interaction (CI) method. One difference from 
the classic CI method is its basis functions, the Slater determinants, being 
over-complete. Another difference is the set of basis states is selected via a 
constrained, importance-sampled random walk. There is not just one set of 
basis functions but many sets. Averaging over these sets is necessary to com- 
pute expectation values. 
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3 Applications 



Except for a very recent application to a nuclear physics model [7], all ap- 
plications of the constrained-path method have been to Hubbard-like models 
of interest to physical chemists and condensed-matter physicists. These mod- 
els can be grouped as one-band, two-band, and three-band models, and each 
group usually has targeted classes of phenomena and materials. The Hub- 
bard models represent considerable simplifications of the complex interactions 
found in actual materials but still represent complicated many-body problems 
for which exact solutions are rare and limited to one- dimensional systems. 
Under these circumstances Monte Carlo methods represent perhaps the only 
controlled means to study the properties of these models and benchmark ap- 
proximate theories for these properties. 

Hubbard models are used to study a sweeping array of intrinsically quantum 
many-electron phenomena. Historically these models were heavily studied for 
possible magnetic phenomena (ferromagnetic, anti-ferromagnetic, paramag- 
netic, etc.). More recently, they have been intensively studied for possible 
representatives of high temperature superconducting materials. In between, 
they were studied as candidate models of heavy fermion and mixed-valence 
materials. 

The principal objects of interest are correlation functions, for example, spin- 
spin, charge- charge, pair-pair correlation functions, because the scaling of 
these functions with systems size gives a measure of symmetry of the ground- 
state. In one and two dimensions, long-range order at finite temperature is gen- 
erally precluded by the Mermin- Wagner theorem. Such states however can ex- 
ist at zero temperature, and their existence is signified by the behavior of cor- 
relation functions with increasing systems size. Long range anti-ferromagnetic 
order for example would be indicated by the (tt, tt) peak of the Fourier trans- 
form of the spin-spin correlation function increasing with increasing system 
size and extrapolating to a non-zero value in the limit of infinite system size. 
Long-range superconducting order would be indicated by the long-range part 
of the pair-pair correlation function remaining a positive number as the sys- 
tems size is extrapolated to infinite volume. Since the sign problem worsens as 
the systems size increases, having a good approximate method, like the CPMC 
method, is an important tool for determining whether long-range order exists 
in these models. 
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3.1 One- Band Models 



The classic one-band Hubbard model is given by the Hamiltonian 

(ij)(T i 

where the double summation is restrict to nearest neighbors, the operators 
and create and destroy an electron of spin a at lattice site i, and — cl^Cia- 
is the number operator at site i. Typically the lattices studied have periodic 
boundary conditions and were square. When the interaction f/ = 0, the model 
is exactly solvable and the non-interacting electrons are described by a single 
band. When the number of electrons equals the number of lattice sites, the 
model is said to be half- filled. At half-fiUing (n) = {Y.a''^ia) = 1- 

In recent years, this model was extensively studied for its magnetic and pos- 
sible superconducting properties. The half-filled model has no sign problem 
for most QMC methods, and QMC studies have played a major role in es- 
tablishing that the ground state of the two-dimensional model has long-range 
anti- ferromagnetic order. This state is consistent with the observed behavior 
of the parent (undoped) state of high temperature superconductors. In these 
materials superconductivity appears when the parent state is doped away from 
half-filling. At dopings relevant to experiment, the sign problem is very bad. 
With the CPMC method it was possible to study the doped Hubbard model 
and to compute various correlation functions as a function of system size. The 
main focus has been on superconducting pairing correlation functions 

To be a bit more definite, the types of correlation functions computed were: 
The spin density structure factor is 

S{k.,,ky) = S{k) = l/NY.e^vi}k-3){sos,), (2) 

3 

where Sj = Uj^ — Uji is the z-component of spin at site j. The charge density 
structure factor is similar to (2), with spin replaced by density, i.e., with the 
— sign in sj replaced by a -|- sign. The electron pairing correlation function is 
defined as 

Pa{j.,Jy)^P{j)^{AUj)A^m, (3) 

where a indicates the symmetry of the pairing. The on-site s-wave pair- 
ing function has As(j) = Cji^Cji, while for rf-wave pairing we used Ad(j) = 
CjiJ2sf{^)(^j+Sij where 6 is (±1,0) and (0, ±1). For 6 along the x-axis, f{6) 
is 1; otherwise it is —1. 
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Fig. 1. Long-range behavior of the (1^2 _y2 pairing correlation function versus distance 
for (n) = 0.85 and a 12 x 12 lattice at [/" = 2, 4, and 8 [8]. This behavior is shown 
for the free-electron electron and CPMC calculations. Also shown is the vertex 
contribution. 

Simulations of the Hubbard model find that the s-wave pairing correlations are 
generally suppressed relative to the d-wave pairing correlations. For a given 
lattice size, increasing the value of U suppresses the d-wave pairing and for 
a given value of U , increasing the lattice size suppresses the d-wave pairing. 
These results are inconsistent with long-range order and are illustrated in 
Figs. 1 and 2 where the d-wave pairing correlation function is shown for a 
12 X 12 and a 16 x 16 lattice as a function of distance In these figures, also 
reported is the "vertex contribution" to the correlation functions defined as 

= P.ij) - PaU) (4) 



where Pa{j) is the contribution of two dressed non- interacting propagators: 
For each term inP^O) of the form (cjcfcjcj), Pa{j) hasaterm like (c|c|) (c|c|). 
When combined. Figs. 1 and 2 indicate the likely absence of long-range super- 
conducting order in the two-dimensional Hubbard model. 

Most projector Monte Carlo ground-state calculations of Hubbard models have 
projected from trial states that were not superconducting. Using a trick, we 
projected from a d-wave BCS superconducting wave function with a BCS 
superconducting order parameter A = 0.5. (Again we use the trial state not 
only for the initial state, but also for importance and constraining states.) 
Within statistical error, we found the same d-wave correlation function as we 
did when we projected from a free-electron wave function. One of our results is 
shown in Fig. 3. This similarity re-enforces the results of Ref. [8] that suggest 
the absence of ODLRO in the two-dimensional Hubbard model. 
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Adding a next-nearest neighbor hopping term of strength t' to the classic 
Hubbard model produces the ff'U model. The CPMC method has been ap- 
plied to this model with the results not being materially different from those 
without the t' . Also adding the t' term, whose presence is suggested by band 
structure calculations, does not seem to change the magnetic properties at 
half-filling, and away from half-filling its addition does not seem to enhance 
superconductivity [9,10]. 

Adding a nearest neighbor interaction, V "n^ii^j with rii — '^io-, produces 
the tUV model. The additional interaction generates additional competing 
effects that leading to a novel co-existence of states [11]. 



3.2 Two-Band Model 



The two-band model is almost always called the periodic Anderson model. In 
this model a d and a / orbital are on each lattice site. These two orbitals per 
unit cell lead to two bands. The Hamiltonian is 

i,a 

where the creation and destruction operators create and destroy d-electrons 
on sites of a square lattice and f-electrons on localized orbitals associated 
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Fig. 2. Long-range behavior of the dx2_y2 pairing correlation function versus distance 
for (n) = 0.85 and a 16 x 16 lattice at {/ = 2 and 4 [8]. This behavior is shown 
for the free-electron electron and CPMC calculations. Also shown is the vertex 
contribution. 
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G O A = Epp^^N =- 1 .09436(4) 

• • A = 0.5 E^p^^ =-1.08975(4) 
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Fig. 3. Pairing correlation functions in the (i^2_y2-wave channel given by the CPMC 
method for a 10 x 10 lattice, U = 4, and (n) = 0.82. The inset shows the long range 
part in detail. The results are the same for the two different trial wave functions: the 

correlations decay quickly with distance [6]. On the other hand the BCS trial wave 
function exhibits ODLRO. Errors bars are smaller than the size of the symbols. 

with these sites. = fj^fi^a is the number operator for f-electrons. Hopping 
only occurs between between neighboring lattice sites {-t term) and between 
a lattice site and its orbital {V term). 

Because of the two bands, the model is half-filled when there are two electrons 
per lattice site. At half-filling the model is said to be symmetric when e/ = 
—U /2. For the symmetric model there is no sign problem and standard QMC 
methods suggest that the model is an insulating anti-ferromagnetic if J7 > 
C/c ~ 2. With the CPMC method studying the magnetic properties of the 
doped model was possible. 

Upon doping with holes, the long-range anti-ferromagnetic order was rapidly 
destroyed. Around three-quarters filling of the lower band, a strong peak devel- 
oped at the (tt, 0) value of the spin structure factor S{{{k) for the /-electrons. 
This peak is shown in Fig. 4, and its development is consistent with the reso- 
nance of two degenerate spin-density waves with wave vectors (tt, 0) and (0, tt) 
[12]. Unestablished is whether this novel state is one of long-range order. 



3.3 Three-Band Model 



The three-band model was constructed with the structure of high-temperature 

superconductors in mind. The common structural feature of these materials 
are Cu02 planes. In addition, these materials possess strong two-dimensional- 
like anisotropy. These properties have focused attention on the physics in these 
planes as being the source of the superconductivity. 
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Fig. 4. Spin-spin correlation function S{[{k) as a function of wavenumber and U for 
|-filling of the lower band [12]. The plot in k-space is along symmetry lines in the 
Brillouin zone 

The three-band model studied by the CPMC method represents the Hamilto- 
nian for the Cu02 planes with only the most relevant Cu and O being kept 



<j,k>a ja j 



(6) 



<i,j> <i,j><y 



In writing the Hamiltonian, we adopted the convention that the operator 
dl,j creates a hole at a Cu 3dx2-y2 orbital and p] creates a hole in an O 
2px or 2py orbital. Ua and Up are the Coulomb repulsions at the Cu and O 
sites respectively, and are the corresponding orbital energies, and Vpd 
is the nearest neighbor Coulomb repulsion. As written, the model has a Cu- 
O hybridization t^j = ^tpd with the minus sign occurring j = i + x/2 and 
j — i — y/2 and also hybridization t^l = ±tpp between oxygen sites with the 
minus sign occurring iox k = j — x/2 — y/2 and k — j + x/2+y /2. The lattice is 
a planar CuO^z structure which has three atoms, one Cu and two O, per unit 
cell. Hence the non-interacting problem has three bands. In a bit of a switch, 
convention for this model has half-filling corresponding to half-filling of the 
anti-bonding (upper) band which is the lowest band in the hole representation. 
When the Hamiltonian is expressed in a hole representation, half-filling then 
corresponds to one hole per unit cell. 

The superconducting pairing correlation functions were studied with the CPMC 
method, and the findings were similar to those found from the studies on the 
Hubbard model: d-wave correlations are stronger than the s-wave ones [13]. 
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Table 1 

Comparison of the exact ground-state energy with the CPMC result with a one 
and two Slater determinant trial wave function for a 2 x 2 system. Parameters are 
e = ep — 64 = 3 and tpp = Up = 0. The use of two Slater determinants as described 
in the text improves the accuracy by an order of magnitude or more [13]. 



Ud CPMC ISD CPMC 2SD Exact 



1 


-5.0613(3) 


-5.0764(2) 


-5.076977 


2 


-4.8475(9) 


-4.8789(7) 


-4.880047 


4 


-4.6073(9) 


-4.6615(6) 


-4.661723 


6 


-4.4884(9) 


-4.5468(6) 


-4.547436 



For a given system size, increasing U suppress the pairing correlations. For a 
given U, increasing the system size suppresses the pairing correlations. 



For the three-band model the binding energy between two holes was also calcu- 
lated. This is a difficult calculation because it requires the accurate calculation 
of the difference between nearly equal energies. It is a significant calculation 
because the binding of holes is a pre-requisite for superconductivity, phase 
separation, or stripe formation. 

To calculate the binding energy for holes, we need to study the half-filled case 
and then the 1 and 2 hole doped cases. In the systems considered (2 x 2, 4 x 2 
and 6x4 unit cells) the 2 hole doped case corresponds to a closed shell case. 
The one hole case is one hole away from a closed shell, and the corresponding 
free-electron wave function is doubly degenerate. In this one-hole case the 
accuracy of the energy is as good as in the closed shell case independently of 
the trial wavefunction used. However, for the half-filled case, which is two holes 
away from a close shell, there are 4 degenerate free-electron states. If we used 
a trial state made from selecting arbitrarily any one of the degenerate states or 
an arbitrary linear combination of these states, the calculated energies would 
not be accurate enough to compute the binding energy. Therefore, we used 
the following procedure for the half-filled case: we diagonalized the interacting 
part of the Hamiltonian in this degenerate subspace, and obtained 2 states 
with energy proportional to Ud and 2 states with zero energy. Of the 2 states 
with zero energy only one of them is a singlet. We used this state which is 
represented by a linear combination of two Slater determinants: (cfcj^|Cfcj^| — 
Ck2,-iCk2,i)\CS) /^/2. In Table 1, we compare the energies obtained using the 
CPMC with the one and two Slater determinants trial wave function and 
energies obtained from exact diagonalization. We see that using two Slater 
determinants improves the accuracy by an order of magnitude or more. The 
accuracy has become better than the closed shell case. 

With this increased accuracy we found parameter ranges where holes bind, 
values of the parameters where binding is optimal, and an increase in binding 
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energy with an increase in system size. Since the appearance of hole binding 
seems decoupled from any enhancement of superconducting correlation, an 
open question is the significance and consequence of this binding. 



4 Parallelization 

If we are considering a system of electrons and lattice sites, the CPU 
time scales as Ny^Nf.N'^ where is the number of walkers. (The number of 
walkers is typically of the order of 200 to 1000 with the larger number usually 
need for the larger values of the interaction strengths.) This number needs to 
be sufficiently large to insure an adequate approximation to the ground state 
and is determined on the basis of experience. The factor NcN'^ comes from 
the scaling of the basic matrix operation that must be performed: The basic 
matrix operations are matrix multiplication, inversion, re-orthogonalization, 
and rank-one updates of matrix inverses. The dominant operation is the matrix 
multiplication propagating a walker. The propagator is represented by a, NxN 
matrix, and each walker by two N x N^j matrices, i.e., one for each spin. 

In general the method is CPU intensive as opposed to memory intensive. As 
a consequence, the basic code usually fits in the memory of one processor and 
the parallelization of the simulation can follow one of several natural paths. 
The simplest path is to give each processor a copy of the code, have it read 
different independent input files, and compute independent runs. A less em- 
barrassing parallelization, and a run time reducing one, is to share the number 
of walkers as equally as possible among as many processors as possible, prop- 
agate the walkers on each processor independently, and combine the results. 
The branched nature of the random walk however requires a slightly different 
procedure. 

Each walker carries a weight, and as the walker propagates, its weight increases 
or decreases. Eventually the large weighted walkers dominate, but propagating 
one of then cost as much computing time as it does to propagate a small 
weighted walker that has little bearing on the final results. Hence carrying the 
small weighted walkers becomes inefficient. In these types of random walks, 
a standard procedure is periodically eliminating walkers with small weights 
and replacing each large weighted walker by many medium weighted walkers 
whose total weight on the average is the same as the single walker. There are 
several schemes to do this, and these procedures are called population control. 
Population control prevents a single walker from ultimately dominating the 
simulation and when used properly reduces the variance of the computed 
results. 

After a population control step, the loads on the processors are unbalanced. 



11 



The natural action then would be to redistribute the load. Unless this load 
represents a relatively large number of walkers per processor, a danger of 
introducing a bias into the simulation exists by performing population control 
on a population too small to be representative. Such a case is expected if 
several hundred walkers distributed across a hundred or so processors. 

Two other options for parallelization are easy to implement. One is to use 
relatively few well-populated processors and have them independently execute 
population control. (This is the procedure we used for small and intermediate 
lattice sizes.) The other, which is more effective in reducing run time for larger 
jobs when a large number of processors are available, is doing population 
control by moving all the walkers onto one processor, performing population 
control there, and then uniformly re-distributing the population. Because the 
amount of information passed between processors is small and the amount 
of time needed for population control is small, this procedure still achieves 
a nearly linear reduction in computation time with the number of processors 
and is simple to implement in any message passing environment. 



5 Concluding Remarks 

The constrained-path method has provided simulators of interacting electron 
systems with a useful tool to study systems sizes impossible by other means. 
With this method it has been possible to investigate the timely and important 
question of the nature of superconducting pairing correlations in candidate 
models for high temperature superconductivity. The results obtained strongly 
suggest the absence of long-range order in these models even though they do 
not provide a rigorous proof of this absence. 

While appearing similar to the fixed-node method long used for continuum 
problem, it is not a fixed node method [14]. Appearing about the same time 
as the CPMC method was a lattice version of the fixed-node method [15,16]. 
The CPMC method appears to give more accurate estimates of quantities like 
the energy and correlation functions, but the fixed-node method can more 
easily propagate some special trial states not expressible as a sum of single 
Slater determinants. 

The concept of constrained random walks has lead to two other quantum 
Monte Carlo methods. One is a constrained random walk at finite tempera- 
tures [10]. The initial results are promising. The significance is that finite tem- 
perature methods have lacked effective analogs even of the fixed-node method 
[17]. The other advance is the constrained phase method. The method al- 
lows for the propagators and Slater determinants to be complex valued: One 
has generalized the sign problem to a phase problem. The ground-state wave 
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functions of a system of electrons in a magnetic field must be complex- valued. 
Making the propagator complex-valued may be convenient if the system has 
long-range forces. Preliminary results on small systems in a magnetic field are 
also promising. 

The 16 X 16 lattices size is most likely the largest for which the simple paral- 
lelization scheme described above is convenient. Often Nc ^ N so the compu- 
tation scales roughly as N^. Thus simulating a 20 x 20 system increases the 
run time by a factor of 4 relative to the 16 x 16 system. This increases the 
simulation time from 2 weeks, for example, to 2 months. It is becoming clear 
that one might need to simulate larger lattice sizes to insure finite-size effects 
are not influencing the results. Going to these larger systems will require dis- 
tributing the matrix operations across many processors. Such operations are 
all fundamental, and procedures for doing this type of parallelization exist. 
We will be exploring their utilization in the near future. 
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